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Abstract 

We investigate the calorimetric liquid-glass transition by performing simulations of a binary 
Lennard-Jones mixture in one through four dimensions. Starting at a high temperature, the 
systems are cooled to T = and heated back to the ergodic liquid state at constant rates. Glass 
transitions are observed in two, three and four dimensions as a hysteresis between the cooling and 
heating curves. This hysteresis appears in the energy and pressure diagrams, and the scanning- 
rate dependence of the area and height of the hysteresis can be described by power laws. The 
one dimensional system does not experience a glass transition but its specific heat curve resembles 
the shape of the D > 2 results in the supercooled liquid regime above the glass transition. As D 
increases, the radial distribution functions reflect reduced geometric constraints. Nearest-neighbor 
distances become smaller with increasing D due to interactions between nearest and next-nearest 
neighbors. Simulation data for the glasses are compared with crystal and melting data obtained 
with a Lennard-Jones system with only one type of particle and we find that with increasing D 
crystallization becomes increasingly more difficult. 
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I. INTRODUCTION 



As a glass-forming liquid is cooled, its relaxation time increases very rapidly and at suf- 
ficiently low temperatures, the relaxation time eventually exceeds the time scale of cooling. 
Thus, provided that crystallization is avoided, the system forms a glass at a temperature T g 



and due to the kinetic effects the non-equili 
history 
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jrium state of this glass depends on its thermal 
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151 ] . In recent years much theoretical effort 



has been made to gain a better understanding of the mechanism(s) responsible for the dra- 
matic slowing down of the relaxation dynamics as well as to investigate the aging dynamics of 



the system once the system has fallen out of equilibrium 16, 



17 



19 



20 



21 



22. 



23|. Most 



of these studies have been done for three dimensional systems since this correspond to the 



ex- 



25]. 



overwhelming majority of real experiments, with some notable exceptions in which the 
perimental) glass transition has been investigated in quasi-two dimensional systems 24j, 
Fewer investigations have been devoted to the question to what extend the phenomenon of 
the glass transition depends on the dimensionality of the system since experimentally it is 
rather difficult to change the dimensions without modifying the interactions. Nevertheless, 
such a study is of interest, since it allows one to estimate the role of the local geometry 
on the glass-forming ability of the system. Although it is evident that with increasing di- 



mensionality D the geometric constraints decrease, it is difficu 
within an analytical calculation in a more quantitative way 26 



t to estimate this tendency 



27 



281 ] . One possibility to 



address the problem is to use computer simulations to determine, with a given interaction 
potential, the dependence of the thermodynamic and structural properties of the system as 
a function of D. 

In the present paper we study therefore how the glass transition of a binary-mixture 
Lennard- Jones (BMLJ) system depends on its dimensionality D. In the past it has been 



shown that such a system shows, for D = 3, many properties of rea 



can serve as a good model for glass-forming systems 10 
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ass formers and thus it 
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38]. 



Here we will focus on the details of the glass transition as characterized by the specific heat 
measurements, we will examine the structure of this system in its glassy state and we will 
compare the results with experimental data. Note that here we will not discuss the relaxation 



dynamics and we refer to Refs. 



39 



40 



41 



44 



451 ] in which such investigations have been 



done for two-dimensional systems. 
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(a) A 80 B 20 (b) A 65 B 35 



FIG. 1: Typical configurations of the A80B20 (a) and A65B35 (b) systems for D = 2 at T = 0. 
These systems have been cooled to T = with cooling rates 7 = —1.0 x 10 -5 and 7 = —1.0 x 10~ 7 , 
(a) and (b) respectively. Open and filled disks represent A and B particles, respectively. Dashed 
lines indicate the (virtual) boundary of the simulation box. 

Section II of this paper discusses the details of the simulation, including the type of system 
and the method by which it is studied. Section III describes the results of the simulations, 
starting with microscopic properties of the system, moving to the macroscopic properties. 
Concluding remarks are given in Section IV. 



II. MODEL AND DETAILS OF THE SIMULATIONS 



Following previous work, we consider a binary mixture of particles A and B, all having 



the same mass, m 



29 



31] . We extend the usual three-dimensional BMLJ molecular 



dynamics simulation to one (D — 1) through four (D = 4) spatial dimensions. A number 
of particles, N, is placed inside a box with edge length L and constant volume V = L D . As 



in previous work, periodic boundary conditions are imposed 46J. The interactions between 
particles are given by the Lennard- Jones potential, U a p(r) = 4e Q/ g[(o" Q/ g/r) 12 — (tfo^/r) 6 ], 
where a,/3 G {A, B }, ctaa = 1.0, €aa = 1-0, oab = 0.8, €ab = 1-5, obb = 0.88 and €bb = 0.5 
29]. Here r is denned through r 2 = Ylk=i x k> where Xk is the k th cartesian component of 



the inter-particle separation. Following common practice, the potential is truncated and 



3 



shifted at r = 2.5a, 



i <7aa being the unit of length, 



291 ]. and Boltzmann's constant, 



a/3 [29|, |46|]. Reduced units are used, wit 
€aa the unit of energy, (ma^ A /48eAA) 1 ^ 2 the unit of time 
ks, is set equal to one. The temperature, T, is controlled by a Nose-Hoover thermostat 46] 
with an effective mass of 48 reduced units. For the molecular dynamics of the particles, the 
equations of motion are integrated using the Verlet algorithm with a time step of 0.02 29]. 
The pressure, P, is monitored using the virial theorem, PV = NT + D" 1 Yli<j % • r^-, where 
fjj is the force and the separation between particles i and j. 

Part of the simulations were carried out for the composition A S pB ? n in three dimensions, 
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37, 



38|. For 



a system that has been studied extensively [29, 
D = 1,2,3, and 4, the composition A 65 B 3 5 was selected because it is, unlike A 80 B2o, stable 
against crystallization for D = 2 at the cooling rates employed here. All simulations begin 
with the system in equilibrium at a sufficiently high initial temperature. The temperature 
is lowered at the rate —7 to T = and then increased back to the initial temperature at 
the rate 7, where 7 = 1.0 x 10~ 3 , 1.0 x 10~ 4 , and 1.0 x 10~ 5 . Further simulations with 
7 = 1.0 x 10 -6 have been performed for D — 2,3, and for D = 2 additional simulations 
were carried out at 7 = 1.0 x 10~ 7 . The A§oB2o system was studied with 7 in the range 
from 7 = 1.0 x 10~ 3 to 7 = 1.0 x 10 -7 . For D = 1,2, and 3 the number of particles 
is A^ = 1000, while A" is 2000 for D — 4. To increase the statistical significance of the 
results, data are averaged over independent runs with different initial configurations. At the 
starting temperatures the relaxation times are very short, equilibrium is rapidly attained 
and the statistically independent starting configurations are readily obtained. The starting 
configurations for D = 1 simulations are random sequences of A and B particles. In general 
the results represent averages over 100 runs. Due to computation time constraints, only 20 
runs are used for the D = 4, A65B35 system as well as for A§oB2o with 7 < 1.0 x 10 -6 . 

The particle density p for the A 80 B 2 o simulations is the same as in previous work with 
L = 9.4 in D = 3, i.e. p = N/L D = 1.204 29(. To establish a common reference point for all 
A65B35 systems at different dimensionalities, the system volume was chosen such that the 
simulation pressure is approximately zero when the temperature reaches T = upon cooling 
at 7 = -1.0 x 10~ 4 . The resulting box edge lengths are L = 1002.5, 29.34, 8.88 (N = 1000) 
and 5.68 (N = 2000) for D = 1,2,3 and 4, respectively. It has been shown that the bulk 



properties of D = 3 BMLJ systems emerge with as few as 65 particles 



43|. 



Further simulations with only A particles (A100) were carried out in order to identify the 



4 



features that distinguish the vitreous from the crystalline state. Box sizes, adjusted such 
that P « for the crystalline state at T = 0, are 1119 (N = 1000), 32.7 (N = 1000), 9.85 
(TV = 1000), and 6.0639 (N = 2048) for D — 1 to 4. Temperature scanning rates ranged 
between 7 = ±1.0 x 10~ 3 and ±1.0 x 10~ 4 . The D = 2 and D = 3 systems crystallize 
spontaneously upon cooling whereas the D = 4 remained in a metastable state. Therefore 
the D = 4 system was heated from an initially prepared fee crystal consisting of 4 4 unit cells 



that each contain 8 particles with 24 equidistant neighbors 



47| 



III. RESULTS AND DISCUSSION 

Since we are interested in the phenomenon of the glass transition, we have to consider the 
stability of the model system against crystallization. Figure [Th. shows a typical configuration 
of the particles for the A 80 B 2 o composition, in two dimensions, cooled to T = at a rate 



of —1.0 x 10~ 5 . As observed previously for a similar system 48[, the configuration contains 
areas of hexagonally crystallized A particles in a matrix of amorphous AB material, i.e. the 
system which is a good glass-former in three dimensions is crystallizing in two dimensions. 
In order to suppress these hexagonal A crystals for D = 2, we selected the composition 
A 65 B 3 5. Figure [Tb shows a typical configuration of a A 65 B3 5 system which has been cooled 
to T = at a rate of —1.0 x 10~ 7 , the slowest rate employed in this study. The structure 
appears to be fully amorphous, as required. In the following we will first discuss the results 
obtained for the A 65 B 3 5 composition. Commonalities and differences between A 65 B 3 5 and 
A 80 B 2 o for D = 3 will be considered at the end of this section. 

In Fig.[2]^a) we show the temperature dependence of the pressure for different values of D. 
The P(T) curves are approximately linear, and their slope increases with the dimensionality. 
Recall that the particle density is adjusted for each dimension such that P ~ as the 
temperature reaches T = upon cooling at 7 = —1.0 x 10~ 4 (see Sec. HT])- Actual pressures 
obtained after cooling to T = at this rate are between —0.10 and 0.00. Cooling the system 
more slowly results in a more relaxed glass state, i.e. a more efficient packing of the particles 
and thus a lower pressure at T = 0. For the rate 7 = —1.0 x 10 -5 (Fig. [2]), pressures range 
from —0.59 (D = 4) to 0.00 (D = 1). Pressures at T = are small compared with pressures 
at T g , with \P(T = 0)| less that 5% of P(T g ), as required for consistency. 

The figure also shows that there is a hysteresis between the cooling curve and the heating 
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FIG. 2: (a) Temperature dependence of the pressure of the Ags!^ systems for D = 1 to D = 4 
upon cooling (dashed lines) followed by heating (solid lines) at 7 = =pl.O x 10 -5 . (b) Temperature 
derivatives of the same pressure data. 

curve for D > 2, which becomes more pronounced as D increases. This hysteresis is seen 
more clearly in the temperature derivative of the pressure, dP/dT (Fig. [2(b)). As it will 
be discussed below, the glass transition temperatures for the data shown in Fig. [2] are 
T g = 0.33, 0.58, and 0.89 for D = 2,3,4, respectively. These values for T„ were obtained 



by first calculating the Active temperature of the system, T/(T), as defined by Tool 49]. 
The procedure is based on analytic approximations for the specific heat of the supercooled 
liquid 50| and the glass (further details are given below). In the liquid state the Active 
temperature equals the temperature, while in the glass state it becomes frozen at a finite 
value. We set T g = T)(0), the limiting value of the Active temperature as the glass stops 
evolving upon cooling to low temperatures. 

Above the temperature range of this hysteresis loop, the system is in equilibrium and 
the pressures upon cooling and heating coincide. A hysteresis is a hallmark of physical 
glasses in the glass transition range [5jJ and the glass transition found here is qualitatively 
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FIG. 3: (a) Radial distribution of A-A pairs for A65B35 cooled to T = at the rate — 1.0 x 10~ 5 . 
Data for D = 1 are multiplied by 0.02. Data for D > 2 are shifted up successively by 5. (b) 
Lennard-Jones potential for A-A interaction, truncated and shifted at r = 2.5. 

consistent with previous observations of a glass transition in cooling curves obtained with a 
D = 3, A80B20 BMLJ system 10J. Last not least we remark that the D — 1 system does not 
show any sign of a glass transition since in one dimension this Lennard-Jones model is not 
sufficiently frustrated to form a glass. (However, other one-dimensional models can show a 



glass transition, see, e.g., Ref. 



52|. 



In order to characterize the local structure of the particles it is useful to consider the 
radial pair distribution functions, g a p(r), defined by [531 ] 



F(jik^£|> (r - r « )) 

ie{a}je{/3} 



where N a is the number of particles of type a, r»j is the distance between particle i and j, 
and the factor F(r,D) normalizes g a p{r) to 1.0 for large r, i.e. F(r, 1) = 2, F(r, 2) = 2nr, 
F(r,3) = 4vrr 2 , and F(r,4) = 2vr 2 r 3 . 

Figures [3^, UK, and[5^ show the A-A, A-B and B-B distributions, respectively, at T = 0. 
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FIG. 4: (a) Radial distribution of A-B pairs for A65B35 cooled to T = at the rate — 1.0 x 10~ 5 . 
Data for D = 1 are multiplied by 0.02. Data for D > 2 are successively shifted up by 10. (b) 
Lennard-Jones potential for A-B interactions, truncated and shifted at r = 2.0. 

For D = 1 one finds a first peak in g a p at a distance a a p2 l / & , i.e. at the location of the 
minimum in the corresponding pair potential (shown in panel b of the figures). This is due 
to the fact that we have adjusted the pressure to be zero at T = 0. The second and third 
A-A peaks in Fig. [3^ for D = 1 are at r = 1.80 and r = 2.25 and correspond to ABA and 
AAA elements of the particle chain, respectively, where the underlined letters indicate the 
atom types considered in the pair distribution. For (?ab(j) for D = 1, chain elements of the 
types ABB and AAB give rise to the peaks at r = 1.89 and r = 2.03, respectively. The 
D = 1 B-B radial distribution peaks are due to, in order of increasing r, BB, BAB and BBB 
chain elements. 

For D = 2, the peak at r = 1.76 corresponds to A-A next-nearest neighbors that are 
separated by a nearest neighbor particle of type B, while the peak at r = 2.20 correspond to 
A particles which are separated by a nearest neighbor particle of type A. The D = 2 peaks 
at r = 1.73 and r = 1.96 in the A-B distribution, Fig. m can be interpreted in a similar 
fashion and correspond to correlations with two A and two B nearest neighbors, respectively. 
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FIG. 5: (a) Radial distribution of B-B pairs for A65B35 cooled to T = at the rate — 1.0 x 10 -5 . 
Data for D = 1 are multiplied by 0.01. Data for D > 2 are successively shifted up by 4. (b) 
Lennard-Jones potential for B-B interactions, truncated and shifted at r = 2.2. 

In <?bb(?") a new B-B peak arises from the long edge of rectangles formed by four B particles 
surrounding an A particle, see Fig. QJb). 

We also mention that a notch can be seen in #ab for D = 2 at r = 2.0, the point where 



46( | . It is unlikely that this 



the force is discontinuous due to the truncation of the potential 
discontinuity has a significant impact on the results. 

The functions g a /3{r) also show a well defined nearest neighbor peak for D = 3 and 4, 
but the peaks at larger distances for these dimensions are much less pronounced than the 
ones found for D = 1 and 2, and it becomes difficult to associate them with particular 
arrangements of particles. Note that also the first nearest neighbor peak becomes broader 
with increasing D since the typical distance between neighboring particles will, at large D, 
be strongly influenced by the type and the number of the particles that are their common 
nearest neighbors. For the A-A correlation the increase of D will, e.g., make it possible 
that two A particles share an increasing number of B particles as first nearest neighbors 
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FIG. 6: Potential energy per particle versus temperature for A65B35 upon cooling (dashed lines) 
followed by heating (solid lines) at 7 = =FlO~ 5 . Potential energies are divided by the number of 
neighbors in a closely packed structure (lignancy = 2, 6, 12 and 24 for D = 1, 2, 3, 4, respectively). 



and, since the A-B interaction is strongly attractive, thus decrease the nearest neighbor 
distance of such an A-A pair. Alternatively an A-A nearest neighbor pair that shares many 
A particles as first nearest neighbors, will have a distance that is somewhat larger than the 
average nearest neighbor A-A distance. 

This mechanism for reducing the geometric constraints between particles that are nearest 
neighbors affects of course also the second, third,... -nearest neighbor configurations. This 
is the reason why with increasing D the radial distribution functions become less structured 
at a given r. In particular the location of the second, third,... -nearest neighbor peaks will 
shift to smaller distances and also the minima between consecutive peaks will be be less 
pronounced. This change in the geometry has the effect that the pressure of the system 
increases since the second nearest neighbor particles move to distances in which the potential 
is steeper/more attractive and hence the virial increases. 

Having discussed the influence of the dimensionality on the structure we now present the 
results regarding the glass transition. One convenient method to investigate this transition 
in real glasses are specific heat measurements. We will consider the temperature dependence 



of the specific heat and energy of the BML J system. (See Refs. jlO, llll, ll2|, ll3|, [jjj, |54|, [55] for 
a discussion how other quantities depend on the cooling rate.) Figure [6] shows the potential 
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FIG. 7: Potential energy per particle versus temperature for systems containing only type A 
particles upon cooling (dashed lines) followed by heating (solid lines) at 7 = =pl0~ 4 for D = 1, 2 
and 3. Data for D = 4 were obtained at 7 = ±10 -3 by heating and melting a closed-packed four- 
dimensional face-centered cubic crystal (lowest curve), followed by cooling and heating. Potential 
energies are divided by the number of neighbors in a closely packed structure (lignancy = 2, 6, 12 
and 24 for D = 1, 2, 3, 4, respectively). 

energy per particle for A 65 B 3 5 as a function of temperature. For ease of comparison the 
potential energies are scaled by the number of nearest neighbors in closed-packed sphere 
structures. These lignancy values are 2, 6, 12 and 24 for D — 1 to D — 4 47J. At T = 0, 
P pa the interaction energy of an isolated pair of A particles is 0.5. From Fig.[6]we recognize 
that for the A65B35 system the actual values are lower due to the higher A-B binding energy 
and the contributions of the next-nearest neighbor interactions. Total energies per particle 
(not shown) are negative, indicating bound states. The slopes of the U (T) /lignancy curves 
decrease with increasing D, reflecting enhanced cohesion and stability of the glass and liquid 
states per neighbor-pair. 

First-order melting and crystallization are distinct from the glass transition. This is 
shown in Fig. [7] for constant-volume systems composed entirely of A particles. Box sizes 
were adjusted such that P ~ for the crystalline state at T = 0, and in order to compare 
the curves for the different values of D we have divided them by the corresponding lignancy. 
(Note that since the binary system is a good glass-former, its crystalline structures are very 
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complex and not really known 56j, |57|].) For D = 1, the T— dependence of the potential 



energy of the chain similar to that of the random A65B35 chain, and no phase transition oc- 
curs for T > 0, as expected. In this case U(T = 0)/2 is 0.5, as expected, since next-nearest 
neighbor interactions are insignificant. The two-dimensional system crystallizes readily upon 
cooling near T = 1.7, and upon reheating melting nearly coincides with the crystallization. 
Upon cooling the D = 3 system stays in a (metastable) supercooled state down to a tem- 
perature T m 0.9 at which it crystallizes. The melting of the resulting crystal is observed 
at around T = 1.8, i.e. at about twice the temperature of crystallization. (Note that the 
exact values for melting and crystallization are probably affected by finite size effects. Fur- 
thermore one should recall that the resulting ordered structure does not have a long range 
positional ordering. However, these two issues are not that relevant here.) 

Interestingly we found that for D = 4 the pure A system would not crystallize sponta- 
neously. Therefore we assembled a dense-packed crystalline face-centered cubic structure 



with 4 4 unit eel 
est neighbors) 



s of eight particles (2048 particles altogether, each with 24 equidistant near- 
. Figure d lowest curve, shows that this structure melts at T m 6.6. 
Once melted, the D = 4 system does not crystallize upon cooling but instead forms a glass 
near T = 1.5 (heating and cooling data are shown). We conclude that with increasing D 
spontaneous crystallization becomes increasingly inhibited, i.e. the glass-forming ability in- 
creases with D. This trend is in agreement with the conclusions we drew in the discussion 
of the radial distribution functions (Figs. [HI HI and H2J), i.e. that the number of possible local 
packings increases rapidly with D and this entropic factor will make crystallization more 
difficult. 

From Figure [7] we also recognize that at T = the potential energy per lignancy decreases 
monotonically with D for liquids as well as for the glasses, in agreement with the trend 
observed for the binary mixture, Fig. [61 This results reflects the fact that the next-nearest 
neighbor atoms move closer to the central atom as D increases and hence lower their potential 
energy. The liquid state is characterized by a strong negative curvature of the U(T) curves, 
while the curvature is much less for the crystal and glass states for D = 2, 3 and 4. 

The specific heat at constant volume, cy = dE/dT, is shown in Fig. [8^. The law of 



Dulong- Petit requires that for classical solids lim cy(T) = D 5§J and we find that the specific 
heat becomes indeed equal to D as temperature approaches zero. Also included in the graph 
is the T— dependence of the specific heat of the A100 crystals. At low T the specific heat of the 
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FIG. 8: (a) Constant volume specific heat, per particle, for A65B35 upon cooling (dashed lines) 
followed by heating (solid lines) at 7 = =Fl.O x 1CP 5 . The dotted lines show the specific heat of 
crystallized A particles upon heating at 7 = +1.0 x 10~ 4 . The cooling data for the crystals (not 
shown) coincide with the data upon heating, (b) Same data as in (a) but now divided by the 
dimension D. 

glasses and and the crystals decreases with increasing temperature. This decrease is due to 
the anharmonicity of the effective potentials near the equilibrium positions of the particles. 
To see this we consider this effect in detail for the D — 1 case. For a chain segment with only 
one type of particle the first three terms of a Taylor expansion of the effective interparticle 
potential are given by U(r) = Uq + sr 2 /2 + Ar 4 , where Uq < is the binding energy at T = 0. 
The quadratic term, with s > 0, leads to the nearly harmonic motion of the particle around 
its equilibrium position that gives rise to the Dulong- Petit result. The leading anharmonic 
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term has a positive coefficient (A > 0) due to the repulsive interaction of the Lennard- Jones 
potential as neighboring particles approach. For a one-particle anharmonic oscillator of mass 
m the classical limit of the quantum-mechanical result 59( is 



C = 1 - 6X(m/s) 2 k B T. (2) 

The negative term shows that the specific heat decreases in the vicinity of T = 0. However, 
the one-particle calculation on which Eq. (T2J) is based overestimates the decrease of the 
specific heat of the D — 1 particle chain by about a factor of two. This discrepancy is due to 
the multi-particle effects in a linear chain. Although it is unfortunately not possible to take 
into account these effects in an exact way, there exist approximation schemes to calculate 
them 
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6l| and it is found that at low T the specific heat does indeed decrease linearly 
with increasing T. These calculations also show that such a T— dependence is only found 
in the specific heat at constant volume, whereas the one at constant pressure increases with 



increasing T 611 ] . 



Figure [H^i shows that for D = 1 the specific heat of the disordered glass coincides within 
the numerical accuracy of the data with the one from the crystal, i.e. the anharmonic effects 
in the two systems are very similar, at least in the T— range considered. The cy(T) data for 
the (D = 1) glass agrees with the one of the corresponding crystal not only in the T— range 
in which there is a linear T— dependence, but also at temperatures at which cy(T) is no 
longer linear. This is thus evidence that these two systems have also similar higher order 
anharmonic effects. 

Similar results are obtained for the case of two dimensions for which we can compare the 
specific heat of the glass at low temperatures with that of a hexagonal crystals consisting 
solely of A particles, Fig. Ek- Type A particles by themselves crystallize readily, typically 
with about 2 vacancies per 1000 A particles, and the specific heat curves upon cooling and 
heating at 7 = =f10 -4 agree with each other (data not shown). We see that up to T = 0.2 
the specific heats of the hexagonal A phase and A65B35 glass nearly coincide. Therefore 
we can conclude that the linear decrease of the specific heat is not just a particularity of 
the glassy state, but instead a general property of both types of condensed Lennard- Jones 
systems at low temperatures. Qualitatively the same results are obtained for D = 3 and 
D = 4. This observation is in agreement with experimental findings since there it has been 
found that it is advantageous to crystallize samples in situ after measuring their specific 
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heats in the supercooled liquid and glassy states. By subtracting the measured specific heat 
of the crystal one obtains the net glass and supercooled liquid signals, and the net glass 
signal is typically indistinguishable from zero 62]. 

The excess specific heat is defined as Acy(T) = Cy qmd (T) — c|! ass (T), where Cy qmd is 
the specific heat of the system on the (metastable) liquid branch and Cy ass is the specific 
heat of the frozen (non-ergodic) glass. At T g we observe a step Acy(T 9 ) as the system 
switches between the two states. This step and the hysteresis between heating and cooling 
characterizes the glass transition. As with the P{T) curves, see Fig. [2j we find no step or 
hysteresis in Acy for D = 1. However, for D > 2 the glass transition can be readily be 
identified, and we see that T g increases as the dimensionality increases. This rise of T g with 
D is consistent with the higher binding energies per nearest neighbor particle pair (Fig. [6]). 
Also the amplitude of the hysteresis increases with D, Fig. [8^,. In Fig. [Hb we normalize 
cy{T) by D, the specific heat at T = 0. Even with this normalization both Acv{T g ) and the 
area of the hysteresis loop increase with D, reflecting the larger number of steric degrees of 
freedom when the particles can move in more dimensions. 

Figure El also shows that above T g (D) the shape of the cy(T) curves are independent of 
D. This similarity of the specific heat suggests that, in terms of their thermal fluctuations, 
D = 1 particle chains behaves like a fluid down to T = 0. This observation is in agreement 
with analytical calculations for soft-sphere systems which find that 2.0 is the minimum 
(fractional) dimension required for a glass transition 261 ] . The absence of a transition for 
D = 1 may be linked to the fact that the A65B35 chains are non-erodic, since the initial 
random order of the particles in the chains remains fixed. Upon cooling these chains cannot 
reach an energetically favorable state with a higher number of nearest A-B neighbors, while 
the systems in higher dimensions can reach such chemically more ordered states. 

Since the glass transition is related to the fact that the system falls out of equilibrium, 
the specific heat curves will, at temperatures around T g , depend on the cooling and heat- 
ing rate. In Fig. [9K we show the specific heat of the A65B35 mixture in three dimensions 
using temperature scanning rates that are varied by three orders of magnitude. As in real 

nnnnnn 

experiments [l|, |2|, |3|, |5|, |7|, |63j, T g decreases with cooling rate. In Fig. [9b we show cy(T) for 



the Ag B 2 o composition, using scanning rates that vary by four decades. Although at a first 
glance the data for the two compositions look quite similar, there are significant differences. 
Firstly in the A65B35 system the difference between the heating and cooling curves is larger 
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FIG. 9: Constant volume specific heat, per particle, of (a) A65B35 and (b) A§oB2o for D = 3 upon 
cooling (dashed lines) followed by heating (solid lines). Curves at rates less than 7 = 10 -3 are 
shown successively shifted up by 0.5. 

than in the case of A 8 oB 2 o- This difference may be related, see the discussion of Fig. [TJ to 
a higher concentration of locally frustrated structures in the A 65 B 3 5 glass, and thus to its 
better glass forming ability and resistance against crystallization. Secondly the A 80 B 2 o data 
has, for the three slowest-rate cooling curves, a unusual feature that is not present in cy(T) 
for the A 65 B35 system in that one sees a small peak in the specific heat that coincides with 
the maximum of the heating curve, see Fig. [9b. This feature could be a sign of incipient 
crystallization, e.g. due to the appearance of sub-critically sized nuclei. 

Since D = 2 simulations require less computational effort that the three dimensional 
systems, slower cooling and heating rates can be investigated. Figure [10] shows the D = 2, 
A65B35 system for scanning rates that vary by four decades and we see that all of them 
have a well defined glass transition. We note a substantial decrease of the width of the 
glass transition, in agreement with the data shown in Fig. [91 In fact, in general there is no 
qualitative difference between this system and the one in D = 3 and hence we can conclude 
that qualitatively the glass transition does not depend on D, if D > 2. 

From Figs. [9] and [10] one recognizes clearly that above T g the specific heat increases 
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FIG. 10: Constant volume specific heat, per particle, for A65B35 in two dimensions upon cooling 
(dashed lines) followed by heating (solid lines). Curves at rates slower than 7 = 10 -3 are shifted 
up successively by 0.2. 

with decreasing temperature, in agreement with previous studies of this Lennard- Jones 

since the relaxation time shows a 



significant non-Arrhenius dependence on T 



29 



30 



3l|, i.e. the system can be considered 



as "fragile" [64j], and experimentally it is known that fragile glass-formers usually have a 



specific heat that increases with decreasing T 65|, |66 |. 



Due to this increase of cy with decreasing T, the excess specific heat Acy(T) 



Cy qmd (T) — Cy ^ s (T) becomes larger at lower temperatures, in accordance with experimen- 



^glass 



tal results [67|. Accordingly the step Acy(T 9 ) becomes larger as T g decreases upon slower 
cooling. Furthermore we see from Figs. [9] and [10] that the heat flow rises above the specific 
heat of the supercooled liquid as the glass regains metastable equilibrium when the system 
is heated to T > T g , an effect that also observed in real experiments and which is due to the 
kinetics of the glass transition. We compare the T— dependence of the specific heat found 
in the present simulations with experimental results by using normalized scales. Following 
Tool, we scale the ordinate by converting the specific heat signal into the Active temperature 
Tf(T) of the system 49j. For this we have linearly extrapolated the specific heat of the glass 
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FIG. 11: Comparison of the D = 3, A65B35, 7 = 10 6 simulation result with the glass transition in 
experimental systems 



51 



68, 



691 ] . For lithium acetate and Ge02 the data are scaled results from 



calorimetric experiments, while for B2O3 and Si02 the curves are based on volume expansion and 
small angle x-ray scattering data, respectively. Temperatures are scaled by T g . Ordinate values 
are scaled by subtracting the (extrapolated) signal of the glass state and dividing the residual by 
the (extrapolated) supercooled liquid signal. Curves other than the BMLJ system are shifted up 
by successive steps of 1. 

to the supercooled liquid regime, and the specific heat of the supercooled liquid was extrap- 
olated into to glass regime 50( . One then subtracts the (extrapolated) glass signal from the 
data as well as from the supercooled liquid curve. Finally one divides the residual curves 
(cooling and heating) by the residual supercooled liquid curve. The scaled data dTf/dT is 
shown as a function of T/T g in Fig. [TTJ Here T g was denned consistently as the value of 
Tf reached upon cooling to the lowest temperature. Note that by construction the scaled 
signal goes to zero at low temperatures and to unity at high temperatures. Also included in 



the graph is the corresponding data from experiments of different glass-formers 



51 



68 



69] 



The comparison of the data from simulation with the one from experimental shows that 
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FIG. 12: The total energy per particle for A65B35 in four dimensions upon cooling (dashed line) 
followed by heating (solid line) at rate 7 = ^1.0 x 10~ 5 . 

the hysteresis at the glass transition of the BMLJ system reproduces indeed the shape ob- 
served in experiments carried out with fragile and strong glass-formers. We also see that the 
width of the glass transition as observed in the simulation is significantly larger than the one 
found even in the strongest glass-forming substances, such as SiC>2, and much larger than 
the fragile glass-former lithium acetate. This difference could have been expected since this 
width depends not only on the fragility of the glass-former but also on the cooling rate 10] , 
and the temperature scanning rates of the simulations are typically 10 10 times faster than 



laboratory rates 



7] 



In the following we consider the dependence of the hysteresis on temperature scanning 
rate, composition, and dimensionality of the system in more detail. Figure [12] shows the 
total energy E(T) upon cooling and heating for the four- dimensional A 65 B 3 5 system. We 
denote the area of the hysteresis loop between heating and cooling curves as A A, and 
their maximum separation as AE. Figure [[3] shows these two quantities as a function 
of temperature scanning rate 7 on logarithmic scales. The data points fall, to a good 
approximation, onto parallel lines, indicating a power-law dependence of AA and AE on 
7. Therefore we fitted the data with the functional form A A (7) = AA (7/10~ 5 ) Q and 
AE(j) = A-E (7/10 -5 ) e , where AA and AE Q are the area and the maximum separation 
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at 7 = 10~ 5 , respectively. The exponents that best fit all data are a = 0.40 ± 0.02 and 
e = 0.17 ± 0.03. Table [J lists the pre-factors AAq and AEq. From this table, and from 
Fig. EU we can conclude that the hysteresis effect increases with increasing D, in agreement 
with the result shown in Fig. [HI i.e. that the glass-forming ability of the system increases 
with its dimensionality. 

A similar analysis can also be done for real systems, although of course only for D = 



68j. It is found that the area and maximum separation of the hysteresis curves also 



show a power-law dependence on the cooling rate, thus showing that the results from the 
present simulations are consistent with real experiments. However, these experiments yield 
exponents of 0.29 for area and 0.15 for the height, i.e. different values from the one found 
here and hence we can conclude that these exponents are not universal, but material-specific 
or characteristic of the regime of temperature scanning rate. 



Since fragile glass- 



brmers systems have relaxation times r(T) that seem to diverge at 



a finite temperature 64j, their effective activation barriers, i.e. the local slope of log(r) 



vs. 1/T is larger than the one for strong glass-formers, if r has macroscopic values, e.g. 
1 second. The temperature range at which the system falls out of equilibrium at the glass- 
transition is therefore trivially related to the fragility of the system, with fragile (strong) 
glass-formers showing a transition in a narrow (wide) temperature range. As mentioned 
above, the present BMLJ system is expected to be a fragile glass-former since its relaxation 
times show a strongly non-Arrhenius T— dependence. Using the power-laws approximations 
for the scanning rate dependence of the hysteresis loop (Fig. IT3l) . we can obtain a rough 
estimate for the behavior of the BMLJ systems at laboratory scanning rates (7 = 10~ 15 ~ 
2K/s, if we identify the A particles as argon atoms). Of course it is highly uncertain whether 
or not the present BMLJ system avoids crystallization at such slow rates. The width of the 
glass transition AT is proportional to the ratio AA/AE and thus its 7-dependence is given 
by the exponent a — e. Referring to Fig. [TT], one can estimate the width of the hysteresis 
loop for A65B35 at 7 = 10 -6 to be AT = 0.5T g . The power-law approximation, with 
a — e ~ 0.23, predicts then for the laboratory rate a AT = 0.004T 9 (or AT = 0.027T g if 
using the experimentally derived values a = 0.29 and e = 0.15). By comparison, the width 
of the (fragile) lithium acetate glass curve in Fig. [11] is about 0.04. We conclude that, based 
on this analysis, the BMLJ system is indeed fragile, in agreement with the data from the 



T— dependence of the relaxation times 29 



30 



311 ] . Last not least we also can conclude 
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System 



AE Q x 10 3 AA n x 10 3 



A65B35, D = 2 



2.1±0.4 



8.4±0.7 



A65B35, D = 3 



6.6±0.8 



31±5 



A 65 B 35 , D = 4 



16.0±1.5 



76±6 



^80620, -D = 3 



4.5±0.2 



18±2 



TABLE I: Values of the pre-factor for the power-law fits AE = A£ (7/10~ 5 ) a17 and AA = 
A J 4o(7/10~ 5 ) ' 40 , representing the height and area of the hysteresis loop in the total energy, re- 
spectively. 

that the fragility depends on the dimensionality of the system, since the width of the glass 
transition increases with increasing D (see Fig. [T3|) . 

IV. SUMMARY 

We have presented the results of molecular dynamics simulations in order to study how 
the calorimetric glass transition of binary Lennard- Jones systems depends on the cooling 
rate and the number of dimensions. The BMLJ systems were cooled at a constant rate —7, 
followed by reheating to the ergodic liquid state at rate 7. We find that the composition 
A-80B2O) a good glass-former in three dimensions, crystallizes in two dimensions if 7 is small, 
whereas the composition A 65 B3 5 does not crystallize for any dimension at the cooling rates 
investigated here. 

For glasses that have been produced with a given (small) 7 we find that the peaks in the 
radial distribution functions become quickly washed out with increasing dimensionality D, 
reflecting fewer geometric constraints. In particular the nearest neighbor peak of the radial 
distribution functions becomes broader, reflecting a wider range of geometric configurations 
in higher dimensions. This is evidence that the glass-forming ability of the system increases 
with increasing dimensionality. As D increases, the first peak in the radial distribution 
function shifts closer to the central atom due to force exerted by the second-nearest neighbors 
on the nearest neighbors, thus resulting in a stronger T— dependence of the pressure. At low 
temperatures the T— dependence of the constant volume specific heat of the glass is very 
close to the one of a one-component Lennard- Jones crystal, showing that the anharmonic 
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rate, y 



FIG. 13: Area, (a), and height, (b), of the hysteresis loop at the glass transition versus the scanning 
rate on logarithmic scales for different compositions and number of spatial dimensions. Symbols 
correspond to the systems presented in the upper left corner of (a). Solid lines are fits to the power 
laws AA = AA)(7/10~ 5 ) a for (a) and AE = A£o(7/10- 5 ) e for (b). 

effects in glasses and crystals are quite similar. 

A glass transition was observed in two, three and four dimensions, whereas no glass 
transition is observed in one dimension. For D = 1 the specific heat curve resembles at all 
temperatures the one of the supercooled liquid for D > 2, and thus we conclude that the 
one dimensional system behaves kinetically like a liquid down to T = 0. 

For the systems that show a glass-transition we find a hysteresis loop (cooling and heating 
cycle) in the energy per particle as well as in the pressure of the system. At a given cooling 
rate the area of this loop and the temperature at which it occurs increase with increasing D. 
Thus this is further evidence that increasing dimensionality raises the glass-forming ability 
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of the system. 

The glass transition becomes sharper with decreasing temperature scanning rate and 
thus the hysteresis loop shrinks. Power-laws can be used to fit scanning rate dependence 
of the area and height of these loops. The exponents describing this 7-dependence seem to 
be independent of composition or dimensionality. A similar analysis of experimental data 
indicate that these exponents are not universal, but appear to be specific for the system 
considered. 
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